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Abstract 

This  paper  addresses  the  modeling  of  nonlinear  constitutive  relations  and  hysteresis  inherent  to  piezo¬ 
ceramic  materials  at  moderate  to  high  drive  levels.  Such  models  are  necessary  to  realize  the  full  potential 
of  the  materials  in  high  performance  control  applications,  and  a  necessary  prerequisite  is  the  develop¬ 
ment  of  techniques  which  permit  control  implementation.  The  approach  employed  here  is  based  on  the 
quantification  of  reversible  and  irreversible  domain  wall  motion  in  response  to  applied  electric  fields.  A 
comparison  with  experimental  data  illustrates  that  because  the  resulting  ODE  model  is  physic-based,  it 
can  be  employed  for  both  characterization  and  prediction  of  polarization  levels  throughout  the  range  of 
actuator  operation.  Finally,  the  ODE  formulation  is  amenable  to  inversion  with  facilitates  the  development 
of  an  inverse  compensator  for  linear  control  design. 


1  Introduction 

Piezoceramic  materials  are  being  employed  in  an  increasing  number  of  structural  and  structural  acoustic 
systems  due  to  the  fact  that  they  are  lightweight,  compact,  relatively  inexpensive  and  provide  the  capability 
for  both  sensing  and  actuating.  These  advantages  are  bolstered  by  the  property  that  they  exhibit  nearly 
linear  dynamics  and  minimal  hysteresis  at  low  drive  levels.  The  restriction  of  the  materials  to  the  quasi- 
linear  range,  however,  limits  their  output  to  levels  which  are  a  fraction  of  those  attainable  at  moderate 
to  high  drive  levels.  The  tradeoff  for  improved  performance  is  the  manifestation  of  the  hysteresis  inherent 
to  the  materials  as  illustrated  in  Figure  1  with  data  from  a  PZT  5A  wafer.  To  utilize  the  full  capability 
of  piezoceramic  actuators,  it  is  necessary  to  model  the  nonlinear  constitutive  relations  and  hysteresis  in  a 
manner  which  is  conducive  to  both  system  design  and  control.  In  this  paper,  we  consider  a  model  based 
on  the  domain  properties  of  the  materials. 

The  model  is  based  on  the  theory  presented  in  [9,  10]  where  it  was  illustrated  for  the  relaxor  ferroelectric 
material  PMN-PT-BT  employed  in  its  ferroelectric  regime.  In  light  of  this  theory,  the  modeling  of  hysteresis 
in  piezoceramic  materials  is  considered  in  two  steps.  In  the  first,  the  hysteresis-free,  or  anhysteretic, 
relation  between  the  applied  field  and  resulting  polarization  is  quantified  through  a  classical  application 
of  Boltzmann  statistics.  This  provides  a  constitutive  relation  which  is  applicable  at  low  drive  levels  when 
hysteresis  is  negligible.  Hysteresis  is  then  modeled  through  the  characterization  of  reversible  and  irreversible 
motion  of  domain  wall  pinned  at  inclusions  in  the  material.  The  combination  of  the  components  provides 
an  ODE  model  which  incorporates  the  nonlinear  constitutive  relations  and  hysteresis  observed  in  various 
PZT  compounds  at  high  drive  levels. 

To  place  this  model  in  the  context  of  existing  hysteresis  models  for  piezoceramic  materials,  we  note 
that  models  can  be  roughly  categorized  as  microscopic,  macroscopic  and  semi-macroscopic  in  nature.  The 
microscopic  theories  focus  on  the  quantification  of  energy  mechanisms  at  the  lattice  or  domain  level.  As  il¬ 
lustrated  in  [7] ,  such  approaches  can  provide  a  high  degree  of  detail  but  are  limited  to  simple  stoichiometries 
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and  are  currently  difficult  to  implement  in  control  design  due  to  the  large  number  of  required  parameters. 
Macroscopic  theories  lie  at  the  opposite  end  of  the  spectrum  and  are  advantageous  when  the  underlying 
physical  mechanisms  are  poorly  understood  or  difficult  to  implement.  This  category  includes  models  based 
on  anhysteretic  curves  shifted  to  coincide  with  remanence  polarizations  [13]  as  well  as  the  Preisach  models 
which  have  been  developed  to  quantify  hysteresis  in  piezoceramic  materials  [1,  2,  3].  While  advantageous  in 
control  design,  these  models  often  require  a  large  number  of  nonphysical  parameters  which  can  limit  their 
utility  in  systems  which  exhibit  transient  dynamics  or  slowly  varying  operating  conditions.  The  model 
employed  here  is  semi-macroscopic  in  the  sense  that  it  is  based  on  energy  relations  for  the  material  but 
employs  macroscopic  averages  to  obtain  a  small  number  (five)  of  effective  parameters  for  certain  material 
properties.  Due  to  the  physical  nature  of  the  parameters,  they  are  easily  obtained  when  specifying  the 
model  for  particular  applications  and  can  be  efficiently  updated  to  accommodate  slowly  changing  operating 
conditions. 

The  model  is  summarized  in  Section  2  and  validated  through  a  comparison  with  data  from  a  PZT  5A 
wafer  in  Section  3.  The  experimental  results  illustrate  both  the  accuracy  of  the  model  and  its  capability 
for  predicting  the  polarization  due  to  changing  drive  levels. 
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Figure  1.  Hysteresis  observed  in  a  PZT  5A  wafer  in  response  to  a  1600  V  input. 


2  Hysteresis  Model 

We  summarize  here  the  theory  presented  in  [9,  10]  for  quantifying  hysteresis  in  ferroelectric  materials. 
This  development  is  considered  in  two  steps.  In  the  first,  the  anhysteretic  polarization  is  modeled  through 
classical  Boltzmann  statistics.  This  provides  the  minimal  energy  state  which  would  be  attained  by  the 
material  in  the  absence  of  pinning  sites  which  restrict  domain  wall  movement.  Hysteresis  is  incorporated 
in  the  second  step  through  the  quantification  of  the  energy  required  to  bend  and  translate  domain  walls. 
This  respectively  provides  reversible  and  irreversible  components  to  the  polarization  which  can  be  summed 
to  obtain  the  total  polarization  for  the  material. 
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2.1  Anhysteretic  Polarization 

The  anhysteretic  polarization  curve  is  that  which  would  be  obtained  in  an  idealized  material  which  is  devoid 
of  inclusions  or  imperfections.  As  illustrated  in  Figure  2,  the  anhysteretic  has  a  burst  region  near  the 
origin,  where  dipole  rotation  toward  preferred  ionic  configurations  produces  large  changes  in  polarization, 
and  exhibits  saturation  at  high  field  levels  when  charge  distributions  prohibit  further  changes.  In  contrast 
to  the  measured  polarization  curves  at  high  drive  levels,  the  anhysteretic  is  single-valued  and  reversible. 

To  quantify  the  dependence  of  the  anhysteretic  polarization  on  field  and  coupling  effects,  we  consider 
the  effective  field  acting  on  moments  in  the  material.  Under  the  assumption  of  fixed  temperatures,  no 
applied  stresses  and  quasi-static  operating  conditions,  the  effective  field  is  taken  to  be 

Ee  =  E  +  aP  (1) 

where  E  and  P  respectively  denote  the  applied  field  and  resulting  polarization,  and  the  term  aP  quantifies 
the  field  due  to  interdomain  coupling.  The  parameter  a  =  Eq/PSi  where  Eq  and  Ps  respectively  denote 
a  scaling  electric  field  and  the  saturation  polarization,  is  typically  estimated  through  a  least  squares  fit 
to  data.  We  note  that  throughout  this  development,  Ps  denotes  the  theoretical  saturation  value  beyond 
which,  polar  interactions  prevent  further  increases  in  polarization.  This  definition  is  in  accordance  with  the 
corresponding  definition  for  the  saturation  magnetization  (e.g.,  see  [5])  and  should  not  be  confused  with 
the  definition  employed  in  several  texts  (e.g.,  [4,  page  36])  for  the  saturation  value  obtained  by  extending 
the  slope  at  tip  reversal  back  through  the  vertical  axis. 

As  detailed  in  [9,  10],  Boltzmann  statistics  is  then  employed  to  quantify  the  probability  of  dipoles 
occupying  certain  energy  states.  The  anhysteretic  models  differ  according  to  the  assumptions  regarding 
the  possible  dipole  orientations  and  dielectric  behavior  of  the  material.  Under  the  assumption  that  the 
material  is  isotropic  and  the  orientation  of  cells  can  be  in  any  direction,  the  balance  of  the  thermal  and 
electrostatic  energies  through  a  Boltzmann  probability  relation  yields  the  Langevin  equation 

Pan  =  Ps  COth  -  —  (2) 

for  the  anhysteretic  polarization.  The  constant  a  is  specified  by  a  =  4^  where  T  and  Tc  are  the  operating 
and  Curie  temperatures  for  the  material.  Because  the  scaling  field  Eq  is  unknown,  this  constant  must  also 
be  estimated  for  a  given  material. 
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Figure  2.  Ising  spin  and  Langevin  models  for  the  anhysteretic  polarization. 
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A  second  model  which  is  commonly  employed  in  ferroelectric  applications  is  the  Ising  spin  relation 

Pan  =  -Pstanh  j  (3) 

which  results  from  the  assumption  that  dipole  moments  can  only  be  oriented  in  the  direction  of  the  field 
or  opposite  to  it.  As  noted  in  Figure  2  where  the  two  models  are  compared,  the  restrictions  on  possible 
moment  orientation  causes  the  Ising  spin  relation  to  saturate  more  quickly  than  the  Langevin  model. 
Further  discussion  concerning  the  merits  of  the  two  relations  as  well  as  details  regarding  their  derivation 
are  provided  in  [9,  10]. 


2.2  Domain  Wall  Model 


The  relations  (2)  or  (3)  can  be  employed  at  low  drive  levels  but  are  inappropriate  at  moderate  to  high 
drive  levels  since  they  do  not  incorporate  the  hysteresis  inherent  to  the  materials.  As  detailed  in  [10]  and 
included  references,  sigmoidal  hysteresis  of  the  type  depicted  in  Figure  1  is  typically  attributed  to  the 
energy  required  to  translate  domain  walls  across  pinning  sites  in  the  material.  The  pinning  of  domain 
walls  by  defects  or  inclusions  in  the  material  is  due  to  the  reduction  in  energy  which  occurs  at  these  sites. 
The  equilibrium  position  for  the  domain  wall  occurs  when  the  potential  energy  is  minimized.  At  low  field 
levels,  the  walls  remain  close  to  the  equilibrium  position  and  the  motion  is  reversible.  From  an  energy 
perspective,  the  variations  are  not  sufficient  to  cross  a  barrier  in  the  potential  well.  The  motion  becomes 
irreversible  when  sufficient  energy  is  provided  to  cross  the  potential  barrier.  Physically,  this  can  occur  when 
the  domain  wall  intersects  a  remote  pinning  site  and  is  the  mechanism  underlying  domain  wall  translations. 
The  resulting  irreversible  polarization  Pirr  and  reversible  polarization  Prev  are  then  summed  to  obtain  the 
total  polarization.  This  approach  follows  that  employed  by  Jiles  and  Atherton  in  their  corresponding 
hysteresis  model  for  ferromagnetic  materials  [6] . 

To  quantify  the  irreversible  polarization,  it  is  noted  in  [9,  10]  that  the  polarization  level  for  a  given 
effective  field  can  be  expressed  as  that  for  the  ideal  case  minus  losses  required  to  break  pinning  sites.  This 
yields  the  relation 

p.  -  p  _  fibdPirr  (A) 

r irr  —  r an  V  v 

where  the  parameter  8  =  sign(dE)  ensures  that  the  energy  required  to  break  pinning  sites  always  opposes 
changes  in  polarization.  The  parameter  k  is  defined  by  k  =  where  n  denotes  the  average  density  of 
pinning  sites,  (S^)  is  the  average  energy  for  180°  walls  and  p  is  an  average  dipole  moment.  Because  the 
density  and  energy  of  individual  pinning  sites  are  unknown,  the  parameter  k  must  be  estimated  for  a  given 
material. 

The  expansion  of  (4)  yields  the  differential  equation 


dPii 


P  —  P 

J  an  irr 


dE 


5k  a(pan  Pirr ) 


specifying  the  irreversible  polarization.  As  discussed  in  [9,  10],  while  this  expressions  is  adequate  in  most 
regimes,  it  can  yield  nonphysical  solutions  when  the  field  is  reversed  from  saturation  for  materials  which 
exhibit  significant  hysteresis  and  are  driven  at  high  levels.  The  enforcement  of  solely  reversible  polarization 
changes  in  this  regime  eliminates  this  discrepancy  and  yields  the  relation 

rip  ~  p  —  p 

irr  ^  1  an  1  irr 

dE  k5  -  a  ( Pan  ~  Pirr ) 

where 

~  [  1  ,  {dE  >  0  and  P  >  Pan}  or  {dE  <  0  and  P  <  Pan} 

5  =  < 

I  0  ,  otherwise . 


4 


The  second  component  of  the  polarization  is  the  reversible  polarization  which  models  the  effects  of 
domain  wall  bending.  To  first  approximation,  this  is  modeled  by  the  relation 

Prev  ~  c{P<m  ~  P%rr )  (6) 

where  c  is  a  parameter  which  must  be  estimated  for  the  specific  application  (see  [9,  10]). 

The  total  polarization  is  then  given  by 


P  —  Prev  +  Pit 


or  equivalently, 

P  =  cPan  +  (1  —  c)Pirr  ■  (?) 

To  implement  the  model,  the  effective  field  for  a  given  field  and  irreversible  polarization  level  is  computed 
using  (1).  This  effective  field  value  is  then  employed  in  either  (2)  or  (3)  to  compute  the  corresponding 
anhysteretic  polarization.  The  subsequent  irreversible  polarization  is  determined  by  numerically  integrating 
(5).  The  total  polarization  is  then  specified  by  (7).  Various  aspects  of  the  model  are  illustrated  in  the  next 
section  through  a  comparison  with  experimental  data. 

3  Model  Validation 

To  illustrate  the  performance  of  the  model,  we  consider  the  characterization  of  the  polarization  generated 
by  a  cylindrical  PZT  5A  wafer.  The  wafer  had  a  diameter  of  1  in  and  was  10  mils  thick.  To  maintain 
quasistatic  operating  conditions,  the  frequency  of  the  input  electric  field  was  taken  to  be  200  mHz.  Because 
the  material  was  initially  polarized  prior  to  the  hysteresis  measurements,  the  first  cycle  was  omitted  because 
it  exhibited  transient  behavior.  Three  complete  steady  state  cycles  were  then  measured  for  input  voltages 
ranging  from  600  V  to  1600  V.  The  corresponding  field  inputs  to  the  model  were  determined  using  the 
relation  E  =  V/d  where  d=  10  mils  is  the  thickness  of  the  wafer. 

Two  strategies  were  considered  to  illustrate  the  capability  of  the  model  to  quantify  the  inherently 
hysteretic  relation  between  the  input  field  and  generated  polarization  at  various  drive  levels.  In  the  first, 
the  parameters  a,  a,  c,  k  and  Ps  were  obtained  through  a  least  squares  fit  to  polarization  data  collected  with 
a  peak  input  voltage  of  600  V.  The  model,  with  these  parameters,  was  then  used  to  predict  the  polarization 
at  drive  levels  up  through  1600  V.  The  second  option  differed  in  that  parameters  were  obtained  by  fitting 
the  1600  V  data  and  the  model  was  then  used  to  predict  the  polarization  behavior  at  lower  drive  levels. 

The  model  behavior  and  corresponding  data  for  the  first  strategy  are  illustrated  in  Figure  3.  The 
parameter  values  a  =  3.7  x  106  Vm/C,  a  =  4.1  x  105  C/m2,  c  =  0.35,  k  =  1.8  x  106  C/m2  and  Ps  =  0.49 
C/m2  were  obtained  through  a  least  squares  fit  to  the  600  V  data  with  the  Langevin  relation  used  to  model 
the  anhysteretic  polarization.  The  model  predictions  at  800  V,  1000  V  and  1600  V  were  then  compared  with 
the  corresponding  data  at  these  drive  levels.  Throughout  the  range  of  drive  levels,  the  model  accurately 
predicts  the  measured  polarization.  We  note  that  the  capability  of  the  model  for  predicting  the  polarization 
at  various  drive  levels  is  due  to  the  fact  that  it  is  based  on  energy  principles. 

The  identification  of  parameters  at  1600  V  produced  model  fits  that  were  slightly  more  accurate  at 
that  level  but  were  less  accurate  when  predicting  at  600  V.  For  this  data  set,  the  600  V  data  appears  to 
provide  more  information  relevant  to  identifying  the  parameters  and  hence  gives  the  model  slightly  better 
predictive  capabilities  than  identification  at  high  drive  levels.  For  general  PZT  applications,  however, 
determination  of  parameters  and  comparison  of  model  performance  at  both  low  and  high  drive  levels  will 
indicate  which  parameters  provide  the  best  overall  fit.  An  advantage  of  this  approach  is  that  once  the 
parameters  are  determined,  the  model  can  be  used  to  quantify  the  polarization  for  a  large  range  of  drive 
levels  simply  by  specifying  the  input  voltage  or  field. 
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4  Concluding  Remarks 


This  paper  summarizes  a  hysteresis  model  for  piezoceramic  actuators  which  is  based  on  domain  wall  theory 
for  ferroelectric  materials  [9,  10].  The  theory  quantifies  the  hysteresis  inherent  to  the  materials  through 
the  consideration  of  the  energy  required  to  bend  and  translate  domain  walls  which  are  pinned  at  defects 
or  inclusions  in  the  material.  This  provides  reversible  and  irreversible  polarization  components  whose  sum 
represents  the  polarization  due  to  an  applied  field.  Characterization  in  this  manner  provides  the  model 
with  the  capability  for  specifying  the  polarization  at  a  variety  of  input  field  levels  with  one  set  of  model 
parameters.  The  flexibility  of  the  model  is  further  augmented  by  the  small  number  (five)  of  required 
parameters  and  the  physical  nature  of  the  parameters  (e.g.,  Ps  is  often  known  a  priori  or  can  be  directly 
obtained  from  the  data). 

The  performance  of  the  model  was  illustrated  through  a  comparison  with  data  from  a  PZT  5A  wafer. 
The  resulting  model  behavior  was  accurate  both  at  the  drive  level  used  to  estimate  the  parameters  as  well 
as  input  fields  up  to  more  than  double  that  magnitude.  This  illustrated  the  flexibility  of  the  model  and  its 
capability  for  prediction  throughout  the  drive  range  of  the  actuator. 

Finally,  the  ODE  nature  of  the  model  makes  it  amenable  to  inversion  through  the  consideration  of 
a  complementary  ODE  in  a  manner  analogous  to  that  described  in  [8].  This  facilitates  the  construction 
of  an  inverse  compensator  which  can  be  used  for  linear  control  design  [11,  12].  The  application  of  these 
techniques  for  linear  control  implementation  for  piezoceramic  actuators  which  exhibit  hysteresis  is  under 
investigation. 
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